Parametric

Assessment of Fit of Transitions

We need to derive appropriate functional forms and define respective survival functions. One reason to favor patametric models is that they can be easier to generalize. Several candidate distributions were considered to model transitions, including Weibull (assume a monotonically increasing or decreasing hazard), Log-logistic (non-monotonic hazards), Gompertz (assume a monotonically increasing or decreasing hazard), Log-normal (non-monotonic hazards), Exponential (constant hazard), Gamma, Generalized gamma & Generalized F (more flexible).

The following plots fitted survival curves from each model (colored lines), with the Kaplan-Meier estimate, in red, obtained from an example of Jackson available here, added to the contributions of Wathers & Cutler available here.


#options(warn=-1)

n_iter<-10000

tiempo_antes_fits<-Sys.time()

#Weathers, Brandon and Cutler, Richard Dr., "Comparison of Survival Curves Between Cox Proportional 
#Hazards, Random Forests, and Conditional Inference Forests in Survival Analysis" (2017). All Graduate 
#Plan B and other Reports. 927. 
#https://digitalcommons.usu.edu/gradreports/927 

#<div style="border: 1px solid #ddd; padding: 5px; overflow-y: scroll; height:350px; overflow-x: scroll; width:100%">            
#https://devinincerti.com/2019/01/01/sim-mstate.html

n_trans <- max(trans_matrix, na.rm = TRUE)
fits_wei <- vector(mode = "list", length = n_trans)
fits_weiph <- vector(mode = "list", length = n_trans)
fits_llogis <- vector(mode = "list", length = n_trans)
fits_gomp <- vector(mode = "list", length = n_trans)
fits_logn <- vector(mode = "list", length = n_trans)
fits_exp <- vector(mode = "list", length = n_trans)
fits_gam <- vector(mode = "list", length = n_trans)
fits_ggam <- vector(mode = "list", length = n_trans)
fits_genf <- vector(mode = "list", length = n_trans)
fits_genf_orig <- vector(mode = "list", length = n_trans)
fits_ggam_orig <- vector(mode = "list", length = n_trans)
fits_rp02 <- vector(mode = "list", length = n_trans)
fits_rp03 <- vector(mode = "list", length = n_trans)
fits_rp04 <- vector(mode = "list", length = n_trans)
fits_rp05 <- vector(mode = "list", length = n_trans)
fits_rp06 <- vector(mode = "list", length = n_trans)
fits_rp07 <- vector(mode = "list", length = n_trans)
fits_rp08 <- vector(mode = "list", length = n_trans)
fits_rp09 <- vector(mode = "list", length = n_trans)
fits_rp10 <- vector(mode = "list", length = n_trans)
fits_c_wei <- vector(mode = "list", length = n_trans)
fits_c_weiph <- vector(mode = "list", length = n_trans)
fits_c_llogis <- vector(mode = "list", length = n_trans)
fits_c_gomp <- vector(mode = "list", length = n_trans)
fits_c_logn <- vector(mode = "list", length = n_trans)
fits_c_exp <- vector(mode = "list", length = n_trans)
fits_c_gam <- vector(mode = "list", length = n_trans)
fits_c_ggam <- vector(mode = "list", length = n_trans)
fits_c_genf <- vector(mode = "list", length = n_trans)
fits_c_genf_orig <- vector(mode = "list", length = n_trans)
fits_c_ggam_orig <- vector(mode = "list", length = n_trans)
fits_c_rp02 <- vector(mode = "list", length = n_trans)
fits_c_rp03 <- vector(mode = "list", length = n_trans)
fits_c_rp04 <- vector(mode = "list", length = n_trans)
fits_c_rp05 <- vector(mode = "list", length = n_trans)
fits_c_rp06 <- vector(mode = "list", length = n_trans)
fits_c_rp07 <- vector(mode = "list", length = n_trans)
fits_c_rp08 <- vector(mode = "list", length = n_trans)
fits_c_rp09 <- vector(mode = "list", length = n_trans)
fits_c_rp10 <- vector(mode = "list", length = n_trans)
km.lc <-vector(mode = "list", length = n_trans)

fits_cox<-list()
fits_c_cox<-list()
#"gengamma" Generalized gamma (stable parameterisation)
#"gengamma.orig"    Generalized gamma (original parameterisation)
#"genf" Generalized F (stable parameterisation)
#"genf.orig"    Generalized F (original parameterisation)
#"weibull"  Weibull
#"gamma"    Gamma
#"exp"  Exponential
#"lnorm"    Log-normal
#"gompertz" Gompertz

library(flexsurv)

#Specify the formula
fitform <- Surv(time, status) ~ 1

for (i in 1:n_trans){  
  fits_wei[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "weibull")
}

for (i in 1:n_trans){  
  fits_weiph[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "weibullph")
}

for (i in 1:n_trans){
  fits_llogis[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "llogis")
}

for (i in 1:n_trans){
  fits_gam[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "gamma")
}
#LOS ERRORES OCURREN CUANDO NO HAY COVARIABLES
#In (function (q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE,  ... :   NaNs produced
#gamma no funcionó: NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced
for (i in 1:n_trans){
  fits_ggam[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "gengamma")
}

for (i in 1:n_trans){
  fits_gomp[[i]] <- flexsurvreg(formula=fitform,
                                data = subset(ms_d_match_surv, trans == i),
                                dist = "gompertz")
}
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109: valor inicial en 'vmmin' no es finito
for (i in 1:n_trans){
  fits_logn[[i]] <- flexsurvreg(formula=fitform,
                                data = subset(ms_d_match_surv, trans == i),
                                dist = "lnorm")
}

for (i in 1:n_trans){
  fits_exp[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "exp")
}

no_attempts <- 20

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
          r <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "gengamma.orig")
      )
    }
    fits_ggam_orig[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
          r <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "genf.orig")
      )
    }
    fits_genf_orig[[i]] <- r
}  
## Warning in flexsurvreg(formula = fitform, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <-  flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "genf")
      )
    }
    fits_genf[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=1,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp02[[i]] <- r
}  
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=2,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp03[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=3,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp04[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=4,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp05[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=5,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp06[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=6,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp07[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=7,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp08[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=8,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp09[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=9,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp10[[i]] <- r
}

for (i in 1:n_trans){
km.lc[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i))
}

transition_label<- plots$title
attr(transition_label,"names")<- plots$trans


get_distinct_hues <- function(ncolor,s=0.5,v=0.95,seed=2125) {
  golden_ratio_conjugate <- 0.618033988749895
  set.seed(seed)
  h <- runif(1)
  H <- vector("numeric",ncolor)
  for(i in seq_len(ncolor)) {
    h <- (h + golden_ratio_conjugate) %% 1
    H[i] <- h
  }
  hsv(H,s=s,v=v)
}
distinct_hues<-get_distinct_hues(20)

layout(matrix(1:n_trans, nc = 2, byrow = F))
for (i in 1:n_trans){
plot(km.lc[[i]], col="red", conf.int=F);
  lines(fits_wei[[i]], col=distinct_hues[1], ci=F);
  lines(fits_weiph[[i]], col=distinct_hues[2], ci=F);
  lines(fits_gomp[[i]], col=distinct_hues[3], ci=F);
  lines(fits_llogis[[i]], col=distinct_hues[4], ci=F);#A0A36D
  lines(fits_gam[[i]], col=distinct_hues[5], ci=F);
  lines(fits_ggam[[i]], col=distinct_hues[6], ci=F);
  lines(fits_ggam_orig[[i]], col=distinct_hues[7], ci=F);
  lines(fits_logn[[i]], col=distinct_hues[8], ci=F);
  lines(fits_exp[[i]],col=distinct_hues[9], ci=F);
  lines(fits_genf[[i]],col=distinct_hues[10], ci=F, lty = 2);
  lines(fits_genf_orig[[i]],col=distinct_hues[11], ci=F, lty = 2);  
  lines(fits_rp02[[i]],col=distinct_hues[12], ci=F, lty = 2);
  lines(fits_rp03[[i]],col=distinct_hues[13], ci=F, lty = 2);
  lines(fits_rp04[[i]],col=distinct_hues[14], ci=F, lty = 2);
  lines(fits_rp05[[i]],col=distinct_hues[15], ci=F, lty = 2);
  lines(fits_rp06[[i]],col=distinct_hues[16], ci=F, lty = 2);
  lines(fits_rp07[[i]],col=distinct_hues[17], ci=F, lty = 3);
  lines(fits_rp08[[i]],col=distinct_hues[18], ci=F, lty = 3);
  lines(fits_rp09[[i]],col=distinct_hues[19], ci=F, lty = 3);
  lines(fits_rp10[[i]],col=distinct_hues[20], ci=F, lty = 3)
  
legend("bottomleft", legend = c("Kaplan-Meier","Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (orig)", "Lognormal", "Exponential", "Generalized F","Generalized F (orig)","RP2","RP3","RP4","RP5", "RP6","RP7", "RP8", "RP9", "RP10"), col = 
         c("red",distinct_hues), 
       title = "Distributions", cex = .95, bty = "n", lty=c(rep(1,10),rep(2,6),rep(3,4)),lwd=3)# lty = 1:2, 
title(main=transition_label[[i]])
}
Figure 19. Vissual Assessment of Survival Curves

Figure 19. Vissual Assessment of Survival Curves

endTime <- Sys.time()

paste0("Time in process: ");endTime - tiempo_antes_fits
## [1] "Time in process: "
## Time difference of 3.497949 mins
#23 min aprox.

#For more complicated models, users should specify what covariate values they want summaries for, rather than relying on the default
#</div>
options(warn=0)

if(no_mostrar==1){
jpeg(paste0(dta_path,"_mult_state_ags/exp_surv_int_only.jpg"), height=30, width= 30, res= 600, units = "in")
layout(matrix(1:n_trans, nc = 2, byrow = F))
for (i in 1:n_trans){
plot(km.lc[[i]], col="red", conf.int=F);
  lines(fits_wei[[i]], col=distinct_hues[1], ci=F);
  lines(fits_weiph[[i]], col=distinct_hues[2], ci=F);
  lines(fits_gomp[[i]], col=distinct_hues[3], ci=F);
  lines(fits_llogis[[i]], col=distinct_hues[4], ci=F);#A0A36D
  lines(fits_gam[[i]], col=distinct_hues[5], ci=F);
  lines(fits_ggam[[i]], col=distinct_hues[6], ci=F);
  lines(fits_ggam_orig[[i]], col=distinct_hues[7], ci=F);
  lines(fits_logn[[i]], col=distinct_hues[8], ci=F);
  lines(fits_exp[[i]],col=distinct_hues[9], ci=F);
  lines(fits_genf[[i]],col=distinct_hues[10], ci=F, lty = 2);
  lines(fits_genf_orig[[i]],col=distinct_hues[11], ci=F, lty = 2);  
  lines(fits_rp02[[i]],col=distinct_hues[12], ci=F, lty = 2);
  lines(fits_rp03[[i]],col=distinct_hues[13], ci=F, lty = 2);
  lines(fits_rp04[[i]],col=distinct_hues[14], ci=F, lty = 2);
  lines(fits_rp05[[i]],col=distinct_hues[15], ci=F, lty = 2);
  lines(fits_rp06[[i]],col=distinct_hues[16], ci=F, lty = 2);
  lines(fits_rp07[[i]],col=distinct_hues[17], ci=F, lty = 3);
  lines(fits_rp08[[i]],col=distinct_hues[18], ci=F, lty = 3);
  lines(fits_rp09[[i]],col=distinct_hues[19], ci=F, lty = 3);
  lines(fits_rp10[[i]],col=distinct_hues[20], ci=F, lty = 3)
  
legend("bottomleft", legend = c("Kaplan-Meier","Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (orig)", "Lognormal", "Exponential", "Generalized F","Generalized F (orig)","RP2","RP3","RP4","RP5", "RP6","RP7", "RP8", "RP9", "RP10"), col = 
         c("red",distinct_hues), 
       title = "Distributions", cex = .95, bty = "n", lty=c(rep(1,10),rep(2,6),rep(3,4)),lwd=3)# lty = 1:2, 
title(main=transition_label[[i]])
}
  dev.off()
}


#
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
fit_flexsurvreg0<-data.frame()
fitted_flexsurvreg0<-data.frame()

dists_no_covs_10s<-cbind.data.frame(covs=c(rep("fits_",(20)*n_trans)),
              formal=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F(original)", paste0("RP0",2:9),"RP10"),1*n_trans),
              dist=c("weibull", "weibullph", "gompertz", "llogis", "gamma", "gengamma","gengamma.orig", "lnorm", "exp", "genf","genf.orig",rep("no dist",9)),
              model=rep(c("wei", "weiph", "gomp", "llogis", "gam","ggam", "ggam_orig", "logn", "exp", "genf","genf_orig", paste0("rp0",2:9),"rp10"),1*n_trans),
              trans=rep(1:n_trans, each=20))

for (i in 1:nrow(dists_no_covs_10s)){  #
  
  cat(paste0("#### Flexible Survival Model (w/ covs): ",
               dists_no_covs_10s[i,"formal"], "; transition: ",dists_no_covs_10s[i,"trans"], "\n \n"))  
    
  model<-paste0("fits_",dists_no_covs_10s[i,"model"])
  
  if(!is.null(get(model)[[dists_no_covs_10s[i,"trans"]]])){  
     fitted_flexsurvreg0<- rbind(fitted_flexsurvreg0,cbind.data.frame(dist=rep(dists_no_covs_10s[i,"formal"],), 
                            trans=rep(dists_no_covs_10s[i,"trans"],),
                            data.table::data.table(summary(get(model)[[dists_no_covs_10s[i,"trans"]]], 
                                                           #newdata= newdat2a, 
                                                           #newtime=unique(fitted_km$time), 
                                                           type = "survival", 
                                                           tidy=T)))) 
    # Generate fit indices
    fit_flexsurvreg0<-rbind(fit_flexsurvreg0,
       cbind(dist= dists_no_covs_10s[i,"formal"],
             transition=dists_no_covs_10s[i,"trans"],
             fitstats.flexsurvreg(get(model)[[dists_no_covs_10s[i,"trans"]]])))
    #the BIC may not be appropriate if none of the candidate models are considered to be close to the ‘true’ model.     
    } else {
    cat(paste0("The model that assumed a ",dists_no_covs_10s[i,"formal"]," distribution for the transition number ",dists_no_covs_10s[i,"trans"]," did not converge \n\n"))
    }
}
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 1
##  
## The model that assumed a Gompertz distribution for the transition number 1 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 2
##  
## The model that assumed a Gompertz distribution for the transition number 2 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 3
##  
## The model that assumed a Gompertz distribution for the transition number 3 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 4
##  
## The model that assumed a Gompertz distribution for the transition number 4 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 4
## 
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
if(no_mostrar==1){
  fit_flexsurvreg0 %>% 
    dplyr::group_by(transition) %>% 
    summarise(mean(ucl,na.rm=T))
  }
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

km.lc0<-list()
fitted_km0<-data.frame()

for (i in 1:n_trans){
  km.lc0[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i))
  fitted_km0<-rbind(fitted_km0,cbind.data.frame(trans=i,fortify(km.lc0[[i]])))
}

#Calculate error
#c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma","Generalized gamma", "Lognormal", "Exponential", "Generalized F")

fitted_flexsurvreg_binned_mix0<-
data.frame(fitted_flexsurvreg0)[,c("dist","trans","time","est","lcl","ucl")] %>% 
  dplyr::left_join(fitted_km0, by=c("trans","time")) %>% 
#there are many observations that was not available due to empirical missing data
#dplyr::filter(!is.na(surv))
  dplyr::group_by(dist,trans) %>% 
  #dplyr::mutate(est= ifelse(is.na(est), mean(est, na.rm=TRUE), est)) %>% 
  #dplyr::mutate(surv= ifelse(is.na(surv), mean(surv, na.rm=TRUE), surv)) %>% 
  tidyr::fill(surv,.direction="updown") %>% 
  tidyr::fill(est,.direction="updown") %>% 
  ungroup()

db_for_apply_rmse0<-
  cbind.data.frame(
      dist=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F (orignial)", paste0("RP0",2:9),"RP10"),n_trans),
      trans=rep(c(1:n_trans),each=20*2))
   
rmse_comp_fits0<- data.frame()
for(i in 1:nrow(db_for_apply_rmse0)){
  rmse<- Metrics::rmse(subset(fitted_flexsurvreg_binned_mix0, dist==db_for_apply_rmse0[i,"dist"] & 
                       trans==db_for_apply_rmse0[i,"trans"])$est,
              subset(fitted_flexsurvreg_binned_mix0, dist==db_for_apply_rmse0[i,"dist"] & 
                       trans==db_for_apply_rmse0[i,"trans"])$surv)

  rmse_comp_fits0<- rbind(rmse_comp_fits0,cbind(dist=db_for_apply_rmse0[i,"dist"],
                                                  residential=db_for_apply_rmse0[i,"residential"],
                                                  trans=db_for_apply_rmse0[i,"trans"],
                                                  rmse=rmse))
}

rmse_comp_fits0_filter<-
  rmse_comp_fits0 %>% 
    dplyr::filter(rmse!="NaN") %>%  
      dplyr::arrange(trans,dist, rmse)%>%
      dplyr::mutate(rmse=round(as.numeric(rmse),4))

setting_type_label<- c(`0`="Outpatient",`1`="Residential") 
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_


We checked the AICs, cumulative relative error (differences between the Kaplan-Meier curve, and the observed behavior of each distribution).


c26 <- c(
  "dodgerblue2", # red
  "green4",
  "#6A3D9A", # purple
  "#FF7F00", # orange
  "gray16", "gold1",
  "skyblue2", "#FB9A99", # lt pink
  "palegreen2",
  "#CAB2D6", # lt purple
  "#FDBF6F", # lt orange
  "gray70", "khaki2",
  "maroon", "orchid1", "deeppink1", "blue1", "steelblue4",
  "darkturquoise", "green1", "yellow4", "yellow3",
  "darkorange4", "brown", "gray40")

fitted_flexsurvreg_binned_mix0_plot<-
fitted_flexsurvreg_binned_mix0 %>% 
  dplyr::mutate(abs((est-surv)/surv)) %>% group_by(dist,trans) %>% dplyr::mutate(cumsum_error=cumsum(`abs((est - surv)/surv)`)) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(text=paste0("Distribution: ",dist,"<br>survival: ",sprintf("%1.2f",est),"<br>Time (in years): ",round(time/365.25,2), "<br>CIs: ",ifelse(!is.na(lcl),"CIs","no CIs"),"<br>Cumsum(abs((est-surv)/surv))= ",round(cumsum_error,2))) %>% 
  dplyr::mutate(text2=paste0("Distribution: KM","<br>survival: ",sprintf("%1.2f",surv),"<br>Time (in years): ",round(time/365.25,2))) %>%
ggplot()+
    geom_line(aes(time, est, color=dist),size=1)+
    geom_line(aes(time, surv), color="red",size=1)+
    scale_color_manual(name="Distributions", values = c26[1:20])+
                         #c("black",brewer.pal(n = 9, name = 'Paired')))+
                         #c("#112A60","#085754","#D3A347","#4F3C91","red","#112A60","#085754","#8F630D","#251363")) +
    facet_wrap(.~trans,labeller = labeller(trans = transition_label), ncol=3, scales = "free_y")+
    theme_bw()+
    theme(legend.position="bottom",
          strip.background = element_rect(fill = "white", colour = "white"))+
  scale_x_continuous(breaks = seq(0, 365.25*12, 2*365.25), labels=seq(0,12,2))+
  #ylim(0,1.25)+
  #theme(axis.text.x = element_blank(), 
  #      panel.grid.major = element_blank(), 
  #      panel.grid.minor = element_blank()) +
  labs(y="",x="")+
  theme(legend.position='none')

m <- list(
  l = 80*1.05,
  r = 80,
  b = 80,
  t = 80*1.05,
  pad = 4
)
p1<-
fitted_flexsurvreg_binned_mix0_plot$data %>% 
  dplyr::mutate(years=round(time/365.25,0)) %>% 
  dplyr::mutate(trans2=trans) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans2=1:n_trans), by="trans2") %>% 
group_by(trans) %>%
group_map(~ plot_ly(data=., x = ~time, y = ~est, color = ~dist, type = "scatter", mode="lines", text=~text) %>% 
            add_trace(y = ~surv, type = 'scatter', mode = 'lines', line = list(color = 'red', width = 1), text=~text2)%>% 
  layout(annotations = list(list(x = 0.5 , y = 1.03, text = ~unique(transition_label), showarrow = F, xref='paper', yref='paper'))) %>% 
  layout(
    xaxis = list(
      ticktext = list(1, 3, 5, 7, 9, 11), 
      textangle = 0,
      tickvals = list(365.25*1, 365.25*3, 365.25*5, 365.25*7, 365.25*9, 365.25*11),
      tickmode = "array"
  )), keep=TRUE) %>%    
subplot(nrows = 3, shareX = T, shareY=T, titleX = F, titleY = F, margin = .05)%>% layout(showlegend = F) %>% #, margin = 0.05) %>% 
  layout(annotations = list(
                list(x = -0.03, y = 0.5, text = "Survival",
                     font = list(color = "black",size = 15),
                     textangle = 270,
                     showarrow = F, xref='paper', yref='paper', size=48,margin =m)))

p1

Figure 20. Vissual Assessment of Survival Curves (w/o covars)

  #htmlwidgets::saveWidget(as_widget(partial_bundle(p1, local=T)), "test1.html")

#https://github.com/ropensci/plotly/issues/1586

                 #%>% 
  #layout(annotations = list(
  #              list(x = 0.5 , y = -0.03, text = "Survival",
  #                   font = list(color = "black",size = 15),
  #                   textangle = 0,
  #                   showarrow = F, xref='paper', yref='paper', size=48,
  #                   automargin = T)
   #             ))
n_trans2 <- max(trans_matrix, na.rm = T)

fit_flexurvreg0_kable<-
fit_flexsurvreg0 %>% 
  dplyr::arrange(transition, AIC) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans_nmb=1:n_trans2),by=c("transition"="trans_nmb")) %>% 
  dplyr::mutate(trans_w_nmb=paste0(sprintf("%02d",transition),")",transition_label)) %>% 
  dplyr::select(trans_w_nmb, dist, Df, n2ll, AIC, AICc, BIC) %>%
  dplyr::mutate_at(vars(n2ll, AIC, AICc, BIC),~ifelse(abs(.)>1e6,NA,format(round(.,2), scientific = FALSE)))
  
 
fit_flexurvreg0_kable %>%   
      knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption = paste0("Table 11. Fit indices of the adjusted survival analyses (intercept-only)"),
               col.names = c("Transition","Distribution", "DF","Negative 2 Log Likelihood","AIC","AICc","BIC"),
               align =c("l",rep('c', 101))) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
  #kableExtra::add_footnote("Note. NA= Null values", notation="none") %>% 
  kableExtra::scroll_box(width = "100%", height = "375px") 
Table 11. Fit indices of the adjusted survival analyses (intercept-only)
Transition Distribution DF Negative 2 Log Likelihood AIC AICc BIC
01)Admission to Readmission RP07 8 120298.90 120314.90 120314.90 120379.05
01)Admission to Readmission RP05 6 120305.11 120317.11 120317.11 120365.22
01)Admission to Readmission RP09 10 120297.95 120317.95 120317.96 120398.14
01)Admission to Readmission RP08 9 120300.57 120318.57 120318.58 120390.74
01)Admission to Readmission RP10 11 120296.87 120318.87 120318.89 120407.08
01)Admission to Readmission Generalized F(original) 4 120312.07 120320.07 120320.07 120352.15
01)Admission to Readmission Generalized F 4 120312.07 120320.07 120320.07 120352.15
01)Admission to Readmission RP06 7 120308.30 120322.30 120322.30 120378.43
01)Admission to Readmission RP04 5 120324.68 120334.68 120334.69 120374.78
01)Admission to Readmission RP03 4 120337.86 120345.86 120345.86 120377.94
01)Admission to Readmission RP02 3 120346.85 120352.85 120352.85 120376.90
01)Admission to Readmission Generalized gamma 3 120557.30 120563.30 120563.31 120587.36
01)Admission to Readmission Lognormal 2 120856.60 120860.60 120860.60 120876.64
01)Admission to Readmission Generalized gamma (original) 3 120982.51 120988.51 120988.52 121012.57
01)Admission to Readmission Log-logistic 2 121344.72 121348.72 121348.73 121364.76
01)Admission to Readmission Weibull (AFT) 2 121686.14 121690.14 121690.14 121706.18
01)Admission to Readmission Weibull (PH) 2 121686.14 121690.14 121690.14 121706.18
01)Admission to Readmission Gamma 2 121798.14 121802.14 121802.14 121818.18
01)Admission to Readmission Exponential 1 122173.19 122175.19 122175.19 122183.21
02)Readmission to Second Readmission Generalized F(original) 4 36233.31 36241.31 36241.32 36268.35
02)Readmission to Second Readmission Generalized F 4 36233.33 36241.33 36241.33 36268.36
02)Readmission to Second Readmission RP05 6 36231.24 36243.24 36243.26 36283.80
02)Readmission to Second Readmission RP06 7 36229.47 36243.47 36243.49 36290.79
02)Readmission to Second Readmission RP04 5 36233.84 36243.84 36243.84 36277.63
02)Readmission to Second Readmission RP08 9 36227.60 36245.60 36245.63 36306.44
02)Readmission to Second Readmission RP07 8 36230.05 36246.05 36246.07 36300.13
02)Readmission to Second Readmission RP09 10 36226.90 36246.90 36246.94 36314.50
02)Readmission to Second Readmission RP10 11 36226.43 36248.43 36248.47 36322.79
02)Readmission to Second Readmission RP03 4 36241.51 36249.51 36249.52 36276.55
02)Readmission to Second Readmission RP02 3 36260.01 36266.01 36266.01 36286.29
02)Readmission to Second Readmission Generalized gamma 3 36330.75 36336.75 36336.75 36357.03
02)Readmission to Second Readmission Lognormal 2 36364.41 36368.41 36368.41 36381.93
02)Readmission to Second Readmission Generalized gamma (original) 3 36395.82 36401.82 36401.82 36422.10
02)Readmission to Second Readmission Log-logistic 2 36482.31 36486.31 36486.31 36499.83
02)Readmission to Second Readmission Weibull (AFT) 2 36604.39 36608.39 36608.39 36621.91
02)Readmission to Second Readmission Weibull (PH) 2 36604.39 36608.39 36608.39 36621.91
02)Readmission to Second Readmission Gamma 2 36630.18 36634.18 36634.18 36647.70
02)Readmission to Second Readmission Exponential 1 36677.51 36679.51 36679.51 36686.27
03)Second to Third Readmission Generalized F(original) 4 NA NA NA NA
03)Second to Third Readmission Generalized F 4 11326.72 11334.72 11334.74 11357.12
03)Second to Third Readmission RP03 4 11328.59 11336.59 11336.61 11358.99
03)Second to Third Readmission RP02 3 11332.32 11338.32 11338.34 11355.12
03)Second to Third Readmission RP04 5 11328.42 11338.42 11338.45 11366.41
03)Second to Third Readmission RP05 6 11327.94 11339.94 11339.98 11373.53
03)Second to Third Readmission RP07 8 11325.20 11341.20 11341.27 11385.99
03)Second to Third Readmission RP08 9 11323.78 11341.78 11341.87 11392.16
03)Second to Third Readmission RP06 7 11327.98 11341.98 11342.04 11381.17
03)Second to Third Readmission RP09 10 11322.00 11342.00 11342.11 11397.99
03)Second to Third Readmission RP10 11 11321.61 11343.61 11343.74 11405.19
03)Second to Third Readmission Generalized gamma 3 11348.80 11354.80 11354.82 11371.60
03)Second to Third Readmission Lognormal 2 11356.81 11360.81 11360.82 11372.01
03)Second to Third Readmission Generalized gamma (original) 3 11367.06 11373.06 11373.07 11389.85
03)Second to Third Readmission Log-logistic 2 11392.05 11396.05 11396.05 11407.24
03)Second to Third Readmission Weibull (AFT) 2 11433.97 11437.97 11437.98 11449.17
03)Second to Third Readmission Weibull (PH) 2 11433.97 11437.97 11437.98 11449.17
03)Second to Third Readmission Gamma 2 11438.21 11442.21 11442.22 11453.41
03)Second to Third Readmission Exponential 1 11440.83 11442.83 11442.84 11448.43
04)Third to Fourth Readmission Generalized F 4 3510.22 3518.22 3518.28 3536.03
04)Third to Fourth Readmission Generalized F(original) 4 3510.22 3518.22 3518.28 3536.03
04)Third to Fourth Readmission RP03 4 3510.40 3518.40 3518.46 3536.21
04)Third to Fourth Readmission RP04 5 3510.24 3520.24 3520.33 3542.50
04)Third to Fourth Readmission RP05 6 3510.10 3522.10 3522.24 3548.83
04)Third to Fourth Readmission RP02 3 3516.30 3522.30 3522.34 3535.66
04)Third to Fourth Readmission Lognormal 2 3519.21 3523.21 3523.23 3532.12
04)Third to Fourth Readmission RP06 7 3509.46 3523.46 3523.63 3554.63
04)Third to Fourth Readmission RP07 8 3508.80 3524.80 3525.03 3560.43
04)Third to Fourth Readmission Generalized gamma 3 3519.20 3525.20 3525.23 3538.56
04)Third to Fourth Readmission RP08 9 3508.39 3526.39 3526.68 3566.47
04)Third to Fourth Readmission Generalized gamma (original) 3 3520.53 3526.53 3526.57 3539.89
04)Third to Fourth Readmission RP09 10 3507.93 3527.93 3528.28 3572.46
04)Third to Fourth Readmission Log-logistic 2 3524.11 3528.11 3528.13 3537.02
04)Third to Fourth Readmission RP10 11 3508.29 3530.29 3530.71 3579.28
04)Third to Fourth Readmission Exponential 1 3536.16 3538.16 3538.17 3542.62
04)Third to Fourth Readmission Weibull (AFT) 2 3535.32 3539.32 3539.34 3548.23
04)Third to Fourth Readmission Weibull (PH) 2 3535.32 3539.32 3539.34 3548.23
04)Third to Fourth Readmission Gamma 2 3535.97 3539.97 3539.99 3548.88


rmse_plot0<-
rmse_comp_fits0_filter %>% 
  dplyr::mutate(trans=as.numeric(trans)) %>% 
  dplyr::mutate(trans2=as.numeric(trans)) %>% 
    dplyr::left_join(cbind.data.frame(transition_label,trans_nmb=1:n_trans2),by=c("trans"="trans_nmb")) %>% 
  dplyr::mutate(trans_w_nmb=paste0(sprintf("%02d",trans),")",transition_label)) %>% 
        
   ggplot()+
  geom_bar(aes(x=dist, y=rmse), position="dodge", stat="identity", alpha=0.4)+
  theme_bw()+
  #geom_errorbar(aes(x=t, ymin=L, ymax=U, color=program), position="dodge", stat="identity", width=0.4,  alpha=0.9, size=1.3)+
  facet_wrap(.~trans_w_nmb, ncol=3, scales="free_y", dir="h") + 
  xlab("") + 
  ylab("")+
  #ylab("State occupancy probabilities") + 
  theme_minimal()+
  theme(legend.position="bottom")+
    theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_blank(),
          legend.box.background = element_blank())+
  theme(axis.text.x = element_text(angle = 90,hjust=0.95,vjust=0.2))#, vjust=0, hjust=0.3

rmse_plot0
Figure 21. Suvival RMSEs, Ten-states Model (w/o covars)

Figure 21. Suvival RMSEs, Ten-states Model (w/o covars)

#install.packages("survHE")
library(survHE)
#First, defines the vector of models to be used
mods <- c("weibull", "weibullph", "llogis", "gamma", "gengamma", "gompertz", "lnorm" ,"exp")#, "genf")

# And then runs the models using MLE via flexsurv
#m2 <- fit.models(formula = formula, data = data, distr = c("exp","wei","lno","llo"), method="inla")
#m3 <- fit.models(formula = formula, data = subset(ms_d_match_surv, trans == 1), distr = mods, method="hmc")
#Error in model.matrix(formula, data)[(mf %>% filter(event == 1))$ID, ] : 
#  subíndice fuera de  los límites

list_aics<-list()
for (i in 1:n_trans){
  m1 <- fit.models(formula = fitform, data = subset(ms_d_match_surv, trans == i), distr = mods)
  model.fit.plot(MLE=m1, stacked = T, #HMC=m3, 
               models = c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Lognormal", "Exponential", "Generalized F"),
               #lab.profile = c("Outpatient","Residential")
               name_legend = "Inferential method")+ 
  ggtitle("Model comparison based on AIC")+
  labs(subtitle=plots[i,"title"])
}


tiempo_antes_fits2<-Sys.time()

newtime0 = seq(from=round(min(ms_d_match_surv$time),0), to= round(max(ms_d_match_surv$time),0))

#_#_#_#_#_#_#_#_#_#_#_#_#_
#covariates

#Specify the formula
fitform2 <- Surv(time, status) ~  factor(tipo_de_plan_res_1)


kernel_haz_est2a<-list()
kernel_haz_est2b<-list()
kernel_haz2a<-list()
kernel_haz2b<-list()
for (i in 1:n_trans){
library("muhaz")
kernel_haz_est2a[[i]] <- muhaz(ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==1),"time"],
                        ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==1),"status"])
kernel_haz2a[[i]] <- data.table(time = kernel_haz_est2a[[i]]$est.grid,
                         est = kernel_haz_est2a[[i]]$haz.est,
                         dist = "Kernel density")
kernel_haz_est2b[[i]] <- muhaz(ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==0),"time"],
                        ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==0),"status"])
kernel_haz2b[[i]] <- data.table(time = kernel_haz_est2b[[i]]$est.grid,
                         est = kernel_haz_est2b[[i]]$haz.est,
                         dist = "Kernel density")
}

haz_int_only2<-
  rbind(cbind(trans=rep(1:n_trans,each=nrow(kernel_haz2a[[i]])),
              tipo_de_plan_res_1=rep(1,nrow(kernel_haz2a[[i]])),
      rbindlist(kernel_haz2a)),
      cbind(trans=rep(1:n_trans,each=nrow(kernel_haz2b[[i]])),
            tipo_de_plan_res_1=rep(0,nrow(kernel_haz2b[[i]])),
      rbindlist(kernel_haz2b)))

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

n_trans2 <- max(trans_matrix, na.rm = T)

dists_w_covs_10s<-cbind.data.frame(covs=c(rep("fits_c_",(20)*n_trans2)),
              formal=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F(original)", paste0("RP0",2:9),"RP10"),1*n_trans2),
              dist=c("weibull", "weibullph", "gompertz", "llogis", "gamma", "gengamma","gengamma.orig", "lnorm", "exp", "genf","genf.orig",rep("no dist",9)),
              model=rep(c("wei", "weiph", "gomp", "llogis", "gam","ggam", "ggam_orig", "logn", "exp", "genf","genf_orig", paste0("rp0",2:9),"rp10"),1*n_trans2),
              trans=rep(1:n_trans2, each=20))

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
no_attempts <- 30

km.lc2a<-list()
km.lc2b<-list()
fits_cox2<-list()
fits_c_cox2<-list()

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "weibull")
    )
  } 
  fits_c_wei[[i]] <- r
}

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "weibullph")
    )
  } 
  fits_c_weiph[[i]] <- r
}

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "llogis")
    )
  } 
  fits_c_llogis[[i]] <- r
}

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "gamma")
    )
  } 
  fits_c_gam[[i]] <- r
}

for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "gengamma")
      )
    }
    fits_c_ggam[[i]] <- r
}

for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i), 
                                   dist = "gompertz")
      )
    }
    fits_c_gomp[[i]] <- r
}  
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109,  : 
##   valor inicial en 'vmmin' no es finito
for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "lnorm")
      )
    }
    fits_c_logn[[i]] <- r
}  

for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "exp")
      )
    }
    fits_c_exp[[i]] <- r
}  


for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "genf")
      )
    }
    fits_c_genf[[i]] <- r
}  


for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
        attempt <- attempt + 1
        try(
            r <- flexsurvreg(formula=fitform2,
                             data = subset(ms_d_match_surv, trans == i),
                             dist = "gengamma.orig")
        )
    }
    fits_c_ggam_orig[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
        attempt <- attempt + 1
        try(
            r <- flexsurvreg(formula=fitform2,
                             data = subset(ms_d_match_surv, trans == i),
                             dist = "genf.orig")
        )
    }
    fits_c_genf_orig[[i]] <- r
}
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Error in optim(method = "BFGS", par = c(mu = 6.94145455507366, sigma = 0.0321172849831113,  : 
##   non-finite finite-difference value [4]
## Warning in flexsurvreg(formula = fitform2, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=1,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp02[[i]] <- r
}  
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=2,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp03[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=3,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp04[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=4,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp05[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=5,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp06[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=6,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp07[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=7,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp08[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=8,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp09[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=9,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp10[[i]] <- r
}

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

fitted_km<-data.frame()
for (i in 1:n_trans2){
km.lc2a[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i & tipo_de_plan_res_1==1))
km.lc2b[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i & tipo_de_plan_res_1==0))

fitted_km<-rbind(fitted_km,cbind.data.frame(trans=i,residential=rep(1,), fortify(km.lc2a[[i]])))
fitted_km<-rbind(fitted_km,cbind.data.frame(trans=i,residential=rep(0,), fortify(km.lc2b[[i]])))
}

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

newdat2a <- data.table::data.table(tipo_de_plan_res_1= factor(c(rep(1,1))))
newdat2b <- data.table::data.table(tipo_de_plan_res_1= factor(c(rep(0,1))))

fitted_flexsurvreg<-data.frame()
fit_flexsurvreg<-data.frame()

for (i in 1:nrow(dists_w_covs_10s)){  #
  
cat(paste0("#### Flexible Survival Model (w/ covs): ",
             dists_w_covs_10s[i,"formal"], "; transition: ",dists_w_covs_10s[i,"trans"], "\n \n"))  
  
model<-paste0("fits_c_",dists_w_covs_10s[i,"model"])

if(!is.null(get(model)[[dists_w_covs_10s[i,"trans"]]])){  
  #Generate databases
 fitted_flexsurvreg<- rbind(fitted_flexsurvreg,cbind.data.frame(dist=rep(dists_w_covs_10s[i,"formal"],), 
                            trans=rep(dists_w_covs_10s[i,"trans"],),
                            residential=rep(1,),
                            data.table::data.table(summary(get(model)[[dists_w_covs_10s[i,"trans"]]], newdata= newdat2a, newtime=unique(fitted_km$time), type = "survival", tidy=T)))) 
  fitted_flexsurvreg<- rbind(fitted_flexsurvreg,cbind.data.frame(dist=rep(dists_w_covs_10s[i,"formal"],), 
                            trans=rep(dists_w_covs_10s[i,"trans"],),
                            residential=rep(0,),
                            data.table::data.table(summary(get(model)[[dists_w_covs_10s[i,"trans"]]], newdata= newdat2b,  newtime=unique(fitted_km$time), type = "survival", tidy=T)))) 
 #t=newtime0, 
 
  # Generate fit indices
  fit_flexsurvreg<-rbind(fit_flexsurvreg,
     cbind(dist= dists_w_covs_10s[i,"formal"],
           transition=dists_w_covs_10s[i,"trans"],
           fitstats.flexsurvreg(get(model)[[dists_w_covs_10s[i,"trans"]]])))
  #the BIC may not be appropriate if none of the candidate models are considered to be close to the ‘true’ model.     
  } else {
  cat(paste0("The model that assumed a ",dists_w_covs_10s[i,"formal"]," distribution for the transition number ",dists_w_covs_10s[i,"trans"]," did not converge \n\n"))
  }
}
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 1
##  
## The model that assumed a Gompertz distribution for the transition number 1 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 1
##  
## The model that assumed a Generalized F(original) distribution for the transition number 1 did not converge 
## 
## #### Flexible Survival Model (w/ covs): RP02; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 4
## 
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
if(no_mostrar==1){
  fit_flexsurvreg %>% 
    dplyr::group_by(trans) %>% 
    summarise(mean(ucl,na.rm=T))
  }

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_


#load("C:/Users/CISS Fondecyt/OneDrive/Escritorio/mult_state_ags.RData")

c26 <- c(
  "dodgerblue2", # red
  "green4",
  "#6A3D9A", # purple
  "#FF7F00", # orange
  "gray16", "gold1",
  "skyblue2", "#FB9A99", # lt pink
  "palegreen2",
  "#CAB2D6", # lt purple
  "#FDBF6F", # lt orange
  "gray70", "khaki2",
  "maroon", "orchid1", "deeppink1", "blue1", "steelblue4",
  "darkturquoise", "green1", "yellow4", "yellow3",
  "darkorange4", "brown", "gray40")

fitted_flexsurvreg_plot<-
fitted_flexsurvreg %>%  
  dplyr::left_join(fitted_km, by=c("time", "trans", "residential")) %>% 
  group_by(trans,dist,residential) %>% 
  tidyr::fill(surv,.direction="down") %>% 
  ungroup() %>% 
  dplyr::mutate(abs((est-surv)/surv)) %>% 
  group_by(trans,dist,residential) %>% 
  dplyr::mutate(cumsum_error=cumsum(replace_na(`abs((est - surv)/surv)`, 0))) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(text=paste0("Distribution: ",dist,"<br>survival: ",sprintf("%1.2f",est),"<br>Time (in years): ",round(time/365.25,2), "<br>CIs: ",ifelse(!is.na(lcl),"CIs","no CIs"),"<br>Cumsum(abs((est-surv)/surv))= ",round(cumsum_error,2),"<br>Residential: ",ifelse(residential==1,"Yes","No"))) %>% 
  dplyr::mutate(text2=paste0("Distribution: KM","<br>survival: ",sprintf("%1.2f",surv),"<br>Time (in years): ",round(time/365.25,2),"<br>Residential: ",ifelse(residential==1,"Yes","No"))) %>%
ggplot()+
    geom_line(aes(time, est, color=dist, linetype=factor(residential)),size=1)+
    geom_line(aes(time, surv, linetype=factor(residential)), color="red",size=1)+
    scale_color_manual(name="Distributions", values = c26[1:20])+
    facet_wrap(.~trans,labeller = labeller(trans = transition_label), ncol=3, scales = "free_y")+
    theme_bw()+
    theme(legend.position="bottom",
          strip.background = element_rect(fill = "white", colour = "white"))+
 # scale_x_continuous(breaks = seq(0, 365.25*12, 2*365.25), labels=seq(0,12,2))+
  labs(y="",x="")+
  theme(legend.position='none')

m <- list(
  l = 80*1.05,
  r = 80,
  b = 80,
  t = 80*1.05,
  pad = 4
)

p2<-
fitted_flexsurvreg_plot$data %>% 
  dplyr::mutate(years=round(time/365.25,0)) %>% 
  dplyr::mutate(trans2=trans) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans2=1:n_trans2), by="trans2") %>% 
group_by(trans) %>%
group_map(~ plot_ly(data=., x = ~time, y = ~est, color = ~dist, type = "scatter", mode="lines", linetype = ~residential, text=~text) %>% 
            add_trace(y = ~surv, type = 'scatter', mode = 'lines', line = list(color = 'red', width = 1, linetype = ~residential), text=~text2)%>% 
  layout(annotations = list(list(x = 0.5 , y = 1.03, text = ~unique(transition_label), showarrow = F, xref='paper', yref='paper'))) %>% 
  layout(
    xaxis = list(
      ticktext = list(1, 3, 5, 7, 9, 11), 
      textangle = 0,
      tickvals = list(365.25*1, 365.25*3, 365.25*5, 365.25*7, 365.25*9, 365.25*11),
      tickmode = "array"
  )), keep=TRUE) %>%    
subplot(nrows = 3, shareX = T, shareY=T, titleX = F, titleY = F, margin = .05)%>% layout(showlegend = F) %>% #, margin = 0.05) %>% 
  layout(annotations = list(
                list(x = -0.1, y = 0.5, text = "Survival",
                     font = list(color = "black",size = 15),
                     textangle = 270,
                     showarrow = F, xref='paper', yref='paper', size=48,margin =m)))
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
#p2

  htmlwidgets::saveWidget(as_widget(p2), "test2_apr22.html")
## Warning in dir.create(target_dir): 'test2_apr22_files\typedarray-0.1' already
## exists
#https://docs.ropensci.org/plotly/reference/partial_bundle.html

File is available in this link.

fit_flexurvreg_kable<-
fit_flexsurvreg %>% 
  dplyr::arrange(transition, AIC) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans_nmb=1:n_trans2),by=c("transition"="trans_nmb")) %>% 
  dplyr::mutate(trans_w_nmb=paste0(sprintf("%02d",transition),")",transition_label)) %>% 
  dplyr::select(trans_w_nmb,dist,Df,n2ll,AIC,AICc,BIC)%>%
  dplyr::mutate_at(vars(n2ll, AIC, AICc, BIC),~ifelse(abs(.)>1e6,NA,format(round(.,2), scientific = FALSE)))  
  
fit_flexurvreg_kable %>%   
      knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption = paste0("Table 12. Fit indices of the adjusted survival analyses"),
               col.names = c("Transition","Distribution", "DF","Negative 2 Log Likelihood","AIC","AICc","BIC"),
               align =c("l",rep('c', 101))) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
  #kableExtra::add_footnote("Note. NA= Null values", notation="none") %>% 
  kableExtra::scroll_box(width = "100%", height = "375px") 
Table 12. Fit indices of the adjusted survival analyses
Transition Distribution DF Negative 2 Log Likelihood AIC AICc BIC
01)Admission to Readmission Generalized F 5 120081.90 120091.90 120091.90 120131.99
01)Admission to Readmission RP07 9 120086.52 120104.52 120104.52 120176.69
01)Admission to Readmission RP05 7 120092.75 120106.75 120106.75 120162.88
01)Admission to Readmission RP09 11 120085.57 120107.57 120107.58 120195.78
01)Admission to Readmission RP08 10 120088.18 120108.18 120108.19 120188.37
01)Admission to Readmission RP10 12 120084.49 120108.49 120108.50 120204.72
01)Admission to Readmission RP06 8 120095.96 120111.96 120111.96 120176.11
01)Admission to Readmission RP04 6 120112.31 120124.31 120124.31 120172.42
01)Admission to Readmission RP03 5 120125.45 120135.45 120135.45 120175.55
01)Admission to Readmission RP02 4 120135.14 120143.14 120143.14 120175.22
01)Admission to Readmission Generalized gamma 4 120283.58 120291.58 120291.58 120323.66
01)Admission to Readmission Lognormal 3 120610.57 120616.57 120616.58 120640.63
01)Admission to Readmission Generalized gamma (original) 4 120750.31 120758.31 120758.31 120790.39
01)Admission to Readmission Log-logistic 3 121114.44 121120.44 121120.45 121144.50
01)Admission to Readmission Weibull (AFT) 3 121486.11 121492.11 121492.11 121516.17
01)Admission to Readmission Weibull (PH) 3 121486.11 121492.11 121492.11 121516.17
01)Admission to Readmission Gamma 3 121605.72 121611.72 121611.72 121635.77
01)Admission to Readmission Exponential 2 121977.73 121981.73 121981.73 121997.76
02)Readmission to Second Readmission Generalized F(original) 5 36233.16 36243.16 36243.17 36276.96
02)Readmission to Second Readmission Generalized F 5 36233.17 36243.17 36243.18 36276.97
02)Readmission to Second Readmission RP05 7 36230.19 36244.19 36244.21 36291.51
02)Readmission to Second Readmission RP06 8 36228.42 36244.42 36244.44 36298.49
02)Readmission to Second Readmission RP04 6 36232.78 36244.78 36244.79 36285.34
02)Readmission to Second Readmission RP08 10 36226.54 36246.54 36246.58 36314.14
02)Readmission to Second Readmission RP07 9 36228.99 36246.99 36247.02 36307.83
02)Readmission to Second Readmission RP09 11 36225.85 36247.85 36247.89 36322.21
02)Readmission to Second Readmission RP10 12 36225.38 36249.38 36249.43 36330.50
02)Readmission to Second Readmission RP03 5 36240.45 36250.45 36250.45 36284.24
02)Readmission to Second Readmission RP02 4 36259.00 36267.00 36267.01 36294.04
02)Readmission to Second Readmission Generalized gamma 4 36330.27 36338.27 36338.27 36365.31
02)Readmission to Second Readmission Lognormal 3 36363.76 36369.76 36369.76 36390.04
02)Readmission to Second Readmission Gompertz 3 36385.89 36391.89 36391.90 36412.17
02)Readmission to Second Readmission Generalized gamma (original) 4 36396.95 36404.95 36404.95 36431.99
02)Readmission to Second Readmission Log-logistic 3 36481.51 36487.51 36487.51 36507.79
02)Readmission to Second Readmission Weibull (AFT) 3 36603.76 36609.76 36609.76 36630.04
02)Readmission to Second Readmission Weibull (PH) 3 36603.76 36609.76 36609.76 36630.04
02)Readmission to Second Readmission Gamma 3 36629.62 36635.62 36635.62 36655.90
02)Readmission to Second Readmission Exponential 2 36677.08 36681.08 36681.08 36694.60
03)Second to Third Readmission Generalized F(original) 5 NA NA NA NA
03)Second to Third Readmission Generalized F 5 11322.92 11332.92 11332.95 11360.91
03)Second to Third Readmission RP03 5 11327.76 11337.76 11337.79 11365.76
03)Second to Third Readmission RP02 4 11331.51 11339.51 11339.53 11361.91
03)Second to Third Readmission RP04 6 11327.60 11339.60 11339.64 11373.19
03)Second to Third Readmission RP05 7 11327.11 11341.11 11341.17 11380.30
03)Second to Third Readmission RP07 9 11324.38 11342.38 11342.47 11392.76
03)Second to Third Readmission RP08 10 11322.94 11342.94 11343.05 11398.92
03)Second to Third Readmission RP06 8 11327.16 11343.16 11343.23 11387.95
03)Second to Third Readmission RP09 11 11321.29 11343.29 11343.42 11404.87
03)Second to Third Readmission RP10 12 11320.82 11344.82 11344.97 11412.00
03)Second to Third Readmission Generalized gamma 4 11345.40 11353.40 11353.42 11375.79
03)Second to Third Readmission Lognormal 3 11354.58 11360.58 11360.59 11377.37
03)Second to Third Readmission Generalized gamma (original) 4 11364.80 11372.80 11372.82 11395.19
03)Second to Third Readmission Gompertz 3 11383.48 11389.48 11389.50 11406.28
03)Second to Third Readmission Log-logistic 3 11390.67 11396.67 11396.68 11413.46
03)Second to Third Readmission Weibull (AFT) 3 11433.30 11439.30 11439.32 11456.10
03)Second to Third Readmission Weibull (PH) 3 11433.30 11439.30 11439.32 11456.10
03)Second to Third Readmission Gamma 3 11437.59 11443.59 11443.60 11460.39
03)Second to Third Readmission Exponential 2 11440.21 11444.21 11444.22 11455.41
04)Third to Fourth Readmission RP03 5 3509.69 3519.69 3519.78 3541.96
04)Third to Fourth Readmission Generalized F(original) 5 3509.95 3519.95 3520.05 3542.22
04)Third to Fourth Readmission Generalized F 5 3509.95 3519.95 3520.05 3542.22
04)Third to Fourth Readmission RP04 6 3509.53 3521.53 3521.66 3548.25
04)Third to Fourth Readmission RP05 7 3509.40 3523.40 3523.58 3554.58
04)Third to Fourth Readmission RP02 4 3515.57 3523.57 3523.64 3541.39
04)Third to Fourth Readmission RP06 8 3508.77 3524.77 3525.00 3560.40
04)Third to Fourth Readmission Lognormal 3 3519.03 3525.03 3525.07 3538.40
04)Third to Fourth Readmission RP07 9 3508.12 3526.12 3526.41 3566.20
04)Third to Fourth Readmission Generalized gamma 4 3519.03 3527.03 3527.09 3544.85
04)Third to Fourth Readmission RP08 10 3507.71 3527.71 3528.06 3572.25
04)Third to Fourth Readmission Generalized gamma (original) 4 3520.26 3528.26 3528.32 3546.07
04)Third to Fourth Readmission RP09 11 3507.24 3529.24 3529.67 3578.23
04)Third to Fourth Readmission Log-logistic 3 3523.61 3529.61 3529.65 3542.97
04)Third to Fourth Readmission Gompertz 3 3524.73 3530.73 3530.77 3544.09
04)Third to Fourth Readmission RP10 12 3507.69 3531.69 3532.20 3585.14
04)Third to Fourth Readmission Exponential 2 3535.32 3539.32 3539.34 3548.23
04)Third to Fourth Readmission Weibull (AFT) 3 3534.49 3540.49 3540.53 3553.86
04)Third to Fourth Readmission Weibull (PH) 3 3534.49 3540.49 3540.53 3553.86
04)Third to Fourth Readmission Gamma 3 3535.13 3541.13 3541.17 3554.49


Session Info

Sys.getenv("R_LIBS_USER")
## [1] "C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2019 (github)/renv/library/R-4.0/x86_64-w64-mingw32;C:/Program Files/R/R-4.0.2/library"
rstudioapi::getSourceEditorContext()$path
## [1] "C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2019 (github)/SUD_CL/Matching_Process15_APR_22.Rmd"
if (grepl("CISS Fondecyt",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_15_apr22.RData")
  } else if (grepl("andre",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("C:/Users/andre/Desktop/SUD_CL/mult_state_15_apr22.RData")
  } else if (grepl("E:",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_15_apr22.RData")
  } else if (grepl("G:",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_15_apr22.RData")
  } else {
    save.image("~/mult_state_15_apr22.RData")
    path.expand("~/mult_state_15_apr22.RData")
  }

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
) %>% 
  DT::datatable(filter = 'top', colnames = c('Row number' =1,'Variable' = 2, 'Percentage'= 3),
              caption = htmltools::tags$caption(
        style = 'caption-side: top; text-align: left;',
        '', htmltools::em('Packages')),
      options=list(
initComplete = htmlwidgets::JS(
      "function(settings, json) {",
      "$(this.api().tables().body()).css({'font-size': '80%'});",
      "}")))